One major task for scientists in his their daily routine is to prepare compelling graphics for technical document, reports or manuscripts for publication. Graphics in form of figures carry the weight of the arguments. They need to be clear, attractive and convincing. according to [@], the difference between good and bad figures always lead into the difference between a highly influential or an obscure paper; a grant/contract won or lost; and a job interview gone well or poorly.

For the last ten years I have squeezed myself into preparing figures for scientific publications and have made throusands of figures. Honestly to say that over this period I have switched from one software to the other in the figure preparation pipeline. I made figures using Microsfoft Excel, OriginPro, SPSS, Matlab, SigmaPlot, matplotlib in python, base R, ggplot2 in R and many others. However, my current preffered tool for making graphics is the ggplot2 package in R. However, looking on the spectrum over the last ten years and the constant switch from one software/tools to the other, I dont expect that I will continue using ggplot2 for the next ten years.

One thing I have learned over the years is that automation should be the main skills for scientists. Automation serve time for data preparation, analysis and producing outputs. In this post I will take you through how to visualize vector field of surface current from drifter observation using the ggplot2 package [@ggplot]. I will further highlight the drawbacks of the default geom of ggplot2 and why it sometimes fail to produce elegant oceanographic plot. Last I will show you how to use altenative geoms from metR package [@metr] to make overcome the challenges inherited in ggplot2 package.

require(metR)
require(tidyverse)
require(lubridate)
require(oce)
require(ocedata)
require(sf)

etopo bathmetry

etopo = sp::read.asciigrid("e:/GIS/ROADMAP/Etopo 1/Tanzania_etopo1/tanz1_-3432.asc")

etopo = etopo %>% as.tibble() %>% rename(lon = 2, lat = 3, bath = 1) %>% filter(bath >= -2400 & bath < 0)


ggplot() +
# metR::geom_contour2(data = etopo, aes(x = lon, y = lat, z = bath),  breaks = seq(-200,0,50), col = "grey70" )+
# metR::geom_text_contour(data = etopo, aes(x = lon, y = lat, z = bath),breaks = seq(-200,0,50), check_overlap = TRUE, rotate = TRUE)+
# metR::geom_contour2(data = etopo, aes(x = lon, y = lat, z = bath),
#                     breaks = seq(0,4000,200))+
  metR::geom_contour_fill(data = etopo, aes(x = lon, y = lat, z = bath),
                          breaks = seq(-4000,0,400), na.fill = TRUE)+
  metR::geom_contour_tanaka(data = etopo, aes(x = lon, y = lat, z = bath),
                            breaks = seq(-4000,0,400))+
  metR::geom_text_contour(data = etopo, aes(x = lon, y = lat, z = bath),
                        breaks = seq(-4000,0,400), check_overlap = TRUE,
                        rotate = TRUE, size = 3)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  labs(x = NULL, y = NULL)

Dataset

The drifter dataset contain surface current information linked to locations in the physical world. This spatial information help us to understand where the high speed current versus low speed current are found in the ocean. It is helpful to visualize this kind of data in the proper geospatial context i.e to show the data on a realistic map. I have filtered the data to cover the western part of the tropical indian ocean. I prepared and arrange drifter information in data frame—a rectangular collection of variables (columns) and observations (rows). The dataset contains observation of surface current worlwide collected by the Global Drifter Program on all major oceans. Among the variables in the dataset are shown in table @ref(tab:tab1) include:

The drifter observations were randomly distributed within the area as shown in figure @ref(fig:fig1) and requires binning—a process of making equal size grid in the area.

ggplot() +
  geom_point(data = drifter.kimbiji %>% sample_frac(1) , aes(x = lon, y = lat))+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = NULL, y = NULL)
The distribution of drifter observation within the area

The distribution of drifter observation within the area

To minimize biasness of sampling, the area was divided into equal size grids showin in figure @ref(fig:fig2).

## convert drifter observation into simple features
drifter.kimbiji.sf = drifter.kimbiji %>% st_as_sf(coords = c("lon", "lat")) %>%
  st_set_crs(4326)

## divide the tropical indian ocean region into grids
drifter.grid = drifter.kimbiji.sf %>% st_make_grid(n = c(50,50))%>%st_sf()

## 
ggplot()+
  geom_sf(data = drifter.grid)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
Gridded area to fill drifter observations

Gridded area to fill drifter observations

Once the area was gridded, then the the mean value of U and V component and the number of observation were calculated in each grid cell.

drifter.kimbiji.sf.se = drifter.kimbiji.sf%>% filter(season=="SE") 

drifter.gridded = drifter.grid %>% 
  mutate(id = 1:n(), contained = lapply(st_contains(st_sf(geometry),drifter.kimbiji.sf.se),identity),
         obs = sapply(contained, length),
         u = sapply(contained, function(x) {median(drifter.kimbiji.sf.se[x,]$u, na.rm = TRUE)}),
         v = sapply(contained, function(x) {median(drifter.kimbiji.sf.se[x,]$v, na.rm = TRUE)}),
         sst = sapply(contained, function(x) {median(drifter.kimbiji.sf.ne[x,]$sst, na.rm = TRUE)})) 

## Then convert the gridded drifter information into data frame

drifter.gridded = drifter.gridded %>% select(obs, u, v,sst) 

## obtain the centroid coordinates from the grid as table
coordinates = drifter.gridded %>% 
  st_centroid() %>% 
  st_coordinates() %>% 
  as_tibble() %>% 
  rename(x = X, y = Y)

## remove the geometry from the simple feature of gridded drifter dataset
st_geometry(drifter.gridded) = NULL

## stitch together the extracted coordinates and drifter information int a single table for SE monsoon season
current.gridded.se = coordinates %>% bind_cols(drifter.gridded) %>% mutate(season = "SE")

drifter.current.gridded = current.gridded.ne %>% bind_rows(current.gridded.se)

After binning, we found that some grids lack the drifter information, therefore, these grids were filled with the U, V values using the Interpolation technique.

## select grids for SE season only
drf.se = drifter.current.gridded %>%select(-sst)%>%filter(season == "SE") %>% na.omit()

## interpolate the U component
u.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$u)
dimension = data.frame(lon = u.se$xg, u.se$zg) %>% dim()

## make a U component data table from interpolated matrix
u.tb = data.frame(lon = u.se$xg, u.se$zg) %>% gather(key = "lata", value = "u", 2:dimension[2]) %>% mutate(lat = rep(u.se$yg, each = dimension[1])) %>% select(lon,lat, u) %>% as.tibble()

## interpolate the V component
v.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$v)

## make the V component data table from interpolated matrix
v.tb = data.frame(lon = v.se$xg, v.se$zg) %>% gather(key = "lata", value = "v", 2:dimension[2]) %>% mutate(lat = rep(v.se$yg, each = dimension[1])) %>% select(lon,lat, v) %>% as.tibble()

## interpolate the observations
ob.se = interpBarnes(x = drf.se$x, y = drf.se$y, z = drf.se$obs)

## make the V component data table from interpolated matrix
obs.tb = data.frame(lon = ob.se$xg, ob.se$zg) %>% gather(key = "lata", value = "obs", 2:dimension[2]) %>% mutate(lat = rep(ob.se$yg, each = dimension[1])) %>% select(lon,lat, obs) %>% as.tibble()

## stitch now the V component intot the U data table and compute the velocity
uv.se = u.tb %>% bind_cols(v.tb %>% select(v), obs.tb %>% select(obs)) %>% mutate(vel = sqrt(u^2+v^2))

Visualising SE monsoon

ggplot() +
  metR::geom_contour_fill(data = uv.se, aes(x = lon, y = lat, z = obs), na.fill = TRUE, bins = 40) + 
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-7,-3.4), xlim = c(38.5, 42))+
  scale_fill_gradientn(name = "Current",colours = oceColorsJet(120), 
                       # limits = c(0,16),
                       # breaks =seq(4,16,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
ggplot() +
  metR::geom_contour_fill(data = uv.se, aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       # limits = c(0,1.6), 
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(max = 1.8, name = expression(speed~(ms^{-1})))
vector field showing the surface current speed and direction

vector field showing the surface current speed and direction 

ggplot()+
metR::geom_streamline(data = uv.se, 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), 
                            alpha = ..step..),
                      L = 2, res = 2, n = 60, 
                      arrow = NULL, lineend = "round")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.5), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = "", y = "")
Massless current flow during the SE season

Massless current flow during the SE season 

## select grids for SE season only
drf.se.sst = drifter.current.gridded %>%filter(season == "SE") %>% na.omit()

## interpolate the U component
sst.se = interpBarnes(x = drf.se.sst$x, y = drf.se.sst$y, z = drf.se.sst$sst)
dimension = data.frame(lon = sst.se$xg, sst.se$zg) %>% dim()

## make sst into a long format data table from interpolated matrix
sst.se.tb = data.frame(lon = sst.se$xg, sst.se$zg) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(sst.se$yg, each = dimension[1])) %>% 
  select(lon,lat, sst) %>% as.tibble()
ggplot() +
  metR::geom_contour_fill(data = sst.se.tb %>% filter(sst > 27 & sst <=30), aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.2), xlim = c(38.7, 40.5))+
  scale_fill_gradientn(name = "Temperature",colours = oceColors9A(120), 
                       # limits = c(0,23),
                       # breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  scale_mag(name = "Speed", max = 1.8, labels = "1.8 m/s")

Visualising NE monsoon

## nelect grids for ne neason only
drf.ne = drifter.current.gridded%>% select(-sst) %>%filter(season == "NE") %>% na.omit()

## interpolate the U component
u.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$u)
dimension = data.frame(lon = u.ne$xg, u.ne$zg) %>% dim()

## make a U component data table from interpolated matrix
u.tb = data.frame(lon = u.ne$xg, u.ne$zg) %>% gather(key = "lata", value = "u", 2:dimension[2]) %>% mutate(lat = rep(u.ne$yg, each = dimension[1])) %>% select(lon,lat, u) %>% as.tibble()

## interpolate the V component
v.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$v)

## make the V component data table from interpolated matrix
v.tb = data.frame(lon = v.ne$xg, v.ne$zg) %>% gather(key = "lata", value = "v", 2:dimension[2]) %>% mutate(lat = rep(v.ne$yg, each = dimension[1])) %>% select(lon,lat, v) %>% as.tibble()

## interpolate the obnervations
ob.ne = interpBarnes(x = drf.ne$x, y = drf.ne$y, z = drf.ne$obs)

## make the V component data table from interpolated matrix
obs.tb = data.frame(lon = ob.ne$xg, ob.ne$zg) %>% gather(key = "lata", value = "obs", 2:dimension[2]) %>% mutate(lat = rep(ob.ne$yg, each = dimension[1])) %>% select(lon,lat, obs) %>% as.tibble()

## stitch now the V component intot the U data table and compute the velocity
uv.ne = u.tb %>% bind_cols(v.tb %>% select(v), obs.tb %>% select(obs)) %>% mutate(vel = sqrt(u^2+v^2))
ggplot() +
  metR::geom_contour_fill(data = uv.ne, aes(x = lon, y = lat, z = obs), na.fill = TRUE, bins = 40) + 
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",colours = oceColorsJet(120), 
                       limits = c(0,23),
                       breaks =seq(4,20,4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))
ggplot() +
  metR::geom_contour_fill(data = uv.ne, aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.ne, aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColorsVelocity(120), 
                       limits = c(0,1.6),
                       breaks =seq(0.2,1.6,.3),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(max = 1.6, name = expression(speed~(ms^{-1})))
Vector field superimposed on surface current velocity showing the speed and direction

Vector field superimposed on surface current velocity showing the speed and direction 

ggplot()+
metR::geom_streamline(data = uv.ne, 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), 
                            alpha = ..step..),
                      L = 2, res = 2, n = 60, 
                      arrow = NULL, lineend = "round")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.6), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+  
  labs(x = "", y = "")
Massless current flow during the NE monsoon

Massless current flow during the NE monsoon 

eacc.a = st_read("eacc_bound.shp")
## Reading layer `eacc_bound' from data source `E:\Data Manipulation\drifter_pemba\eacc_bound.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 39.2514 ymin: -6.619414 xmax: 40.16052 ymax: -4.362925
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
drifter.eacc = drifter.kimbiji.sf %>% st_intersection(eacc.a)
st_geometry(drifter.eacc) = NULL

drifter.monthly = drifter.eacc %>%
filter( year >=1994) %>%
group_by(month, year) %>%
summarise(u = median(u, na.rm = T),
v = median(v, na.rm = T),
vel = sqrt(u^2+v^2))
ggplot()+
  metR::geom_contour_fill(data = drifter.monthly,
  aes(x = year, y = month, z = vel), na.fill = TRUE, bins = 70)+
  metR::geom_contour2(data = drifter.monthly,
  aes(x = year, y = month, z = vel), na.fill = TRUE)+
  metR::geom_text_contour(data = drifter.monthly,
  aes(x = year, y = month, z = vel), na.fill = TRUE, check_overlap = TRUE, skip = 2)+
  scale_fill_gradientn(colours = oce.colorsVelocity(120), breaks = seq(0.2,1.6,.3))+
  scale_y_reverse(breaks = 1:12, 
                  label = month(seq(dmy(010119),  dmy(311219), by = "month"), 
                                label = TRUE, abbr = TRUE)) +
  scale_x_continuous(breaks = seq(1997,2018,4))+
  # theme_void() +
  theme(legend.key.height = unit(1.4, "cm"),
        axis.text = element_text(size = 12, colour = 1))+
  labs(x = NULL, y = NULL)

sst.anomaly = current.gridded.se %>%select(-season) %>% 
  bind_cols(current.gridded.ne %>% select(sst.ne = sst)) %>% 
  mutate(sst.diff = sst.ne - sst) %>% 
  filter(!is.na(sst.diff))

# interpolate the U component
sst.an = interpBarnes(x = sst.anomaly$x, y = sst.anomaly$y, z = sst.anomaly$sst.diff)
dimension = data.frame(lon = sst.an$xg, sst.an$zg) %>% dim()

## make sst into a long format data table from interpolated matrix
sst.an.tb = data.frame(lon = sst.an$xg, sst.an$zg) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(sst.an$yg, each = dimension[1])) %>% 
  select(lon,lat, sst) %>% as.tibble()

uv.an = uv.ne %>%
  bind_cols(uv.se %>% select(u.se = u, v.se = v)) %>% 
  mutate(u.an = u.se-u, v.an = v.se-v)


ggplot() +
  metR::geom_contour_fill(data = sst.an.tb, 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.an, aes(x = lon, y = lat, dx = u.an, dy = v.an), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .75, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.2), xlim = c(38.7, 40.5))+
  # scale_fill_gradientn(name = "Temperature",colours = oceColors9A(120), 
  #                      # limits = c(0,23),
  #                      # breaks =seq(4,20,4),
  #                      na.value = "white")+
  metR::scale_fill_divergent(midpoint = 0)+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  scale_mag(name = "Speed", max = 1.2, labels = "1.8 m/s")+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

vector.se =ggplot() +
  metR::geom_contour_fill(data = uv.se %>% filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2), 
                          aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.se%>%filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2),  
                    aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.5), xlim = c(39.25, 40.2))+
  scale_fill_gradientn(name = expression(ms^{-1}),
                       colours = oceColorsVelocity(120), 
                       limits = c(0,1.8),
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.1, "cm"), 
        legend.background = element_blank(),
        axis.text.y = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(name = "", max = 1.8, labels = "1.8 m/s")+
  scale_x_continuous(breaks = c(39.4,40.1))+
  scale_y_continuous(breaks = c(-5.8, -4.6))


vector.ne = ggplot() +
  metR::geom_contour_fill(data = uv.ne%>%filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2), 
                          aes(x = lon, y = lat, z = vel), na.fill = TRUE, bins = 70) + 
  metR::geom_vector(data = uv.ne%>%filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2),  
                    aes(x = lon, y = lat, dx = u, dy = v), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.5), xlim = c(39.25, 40.2))+
  scale_fill_gradientn(name = expression(speed~(ms^{-1})),
                       colours = oceColorsVelocity(120), 
                       limits = c(0,1.8),
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag(max = 1.8, name = expression(ms^{-1}))+
  scale_x_continuous(breaks = c(39.4,40.1))+
  scale_y_continuous(breaks = c(-5.8, -4.6))

egg::ggarrange(vector.ne, vector.se, ncol = 2)

streamline.ne = ggplot()+
metR::geom_streamline(data = uv.ne %>%filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2), 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), alpha = ..step..),
                      #arrow = NULL, 
                      #lineend = "round"
                      L = .75, res = 2, n = 40)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.5), xlim = c(39.25, 40.2))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.5), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+
  scale_x_continuous(breaks = c(39.4,40.1))+
  scale_y_continuous(breaks = c(-5.8, -4.6))+  
  labs(x = "", y = "")

streamline.se = ggplot()+
metR::geom_streamline(data = uv.se%>%filter(lon >= 39.2 & lon < 40.3 & lat > -6.4 & lat < -4.2), 
                        aes(x = lon, y = lat, dx = u, dy = v, 
                            color = sqrt(..dx..^2 + ..dy..^2), alpha = ..step..),
                      #arrow = NULL, 
                      #lineend = "round"
                      L = .75, res = 2, n = 40)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.0,-4.5), xlim = c(39.25, 40.2))+
  scale_color_viridis_c(guide = "none")+
  scale_size(range = c(0, 1.5), guide = "none") +
  scale_alpha(guide = "none") +
  theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1),
        axis.text.y = element_blank())+
  scale_x_continuous(breaks = c(39.4,40.1))+
  scale_y_continuous(breaks = c(-5.8, -4.6))+  
  labs(x = "", y = "")

egg::ggarrange(streamline.ne, streamline.se, ncol = 2)

Wind

load("e:/Data Manipulation/xtractomatic/quikscat_wio_2006.RData")


wind.x = wind_x[["data"]]
wind.y = wind_y[["data"]]
date = wind_x$time %>% as.Date()
lon = wind_x$longitude
lat = wind_x$latitude
wind.x.tb = list()
wind.y.tb = list()

for (j in 1:length(date)){
wind.x.tb[[j]] =  data.frame(lon, wind.x[,,j] %>% as.data.frame()) %>%
  gather(key = "lata", value = "wind_x", 2:362) %>%
  mutate(lat = rep(lat, each = 321), time = date[j]) %>%
  filter(lon >=38 & lon <= 42 & lat >=-8 & lat <= -3)

wind.y.tb[[j]] =  data.frame(lon, wind.y[,,j] %>% as.data.frame()) %>%
  gather(key = "lata", value = "wind_y", 2:362) %>%
  mutate(lat = rep(lat, each = 321), time = date[j])%>%
  filter(lon >=38 & lon <= 42 & lat >=-8 & lat <= -3)
}
wind.x.tb = wind.x.tb %>% bind_rows()
wind.y.tb = wind.y.tb %>% bind_rows()

wind.data = wind.x.tb %>% bind_cols(wind.y.tb %>% select(wind_y))
wind.data = wind.x.tb %>% 
  bind_cols(wind.y.tb)  %>%
  mutate(month = month(time),
         season = month,
         season = replace(season,season %in% c(11,12,1,2,3, 4), "NE"),
         # season = replace(season,season %in% c(4,10), "IN"),
         season = replace(season,season %in% c(5,6,7,8,9,10 ), "SE"))%>%
  select(lon, lat, time,month, season, wind_x, wind_y)
wind.data.season = wind.data %>% 
  group_by(lon,lat, season) %>%
  summarise(wind.x = mean(wind_x, na.rm = TRUE),
            wind.y = mean(wind_y, na.rm = TRUE),
            speed = sqrt(wind.x^2 + wind.y^2))

wind.data.monthly = wind.data %>% 
  group_by(lon,lat, month) %>%
  summarise(wind.x = mean(wind_x, na.rm = TRUE),
            wind.y = mean(wind_y, na.rm = TRUE),
            speed = sqrt(wind.x^2 + wind.y^2))
wind.se = ggplot()+
  metR::geom_vector(data = wind.data.monthly%>%filter(month == 8), 
                    aes(x = lon, y = lat, dx = wind.x, dy = wind.y), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .6, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.height = unit(.25, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1), 
        axis.text.y = element_blank(),
        legend.key = element_blank())+ 
  scale_mag(name = expression(speed~(ms^{-1})), max = 15)+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))



wind.ne = ggplot() +
  metR::geom_vector(data = wind.data.monthly%>%filter(month == 1), 
                    aes(x = lon, y = lat, dx = wind.x, dy = wind.y), 
                    arrow.angle = 20, arrow.type = "open", arrow.length = .6, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = c(.15,.85),
        legend.key.height = unit(.25, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1),
        legend.key = element_blank())+ 
  scale_mag(name = expression(ms^{-1}), max = 15)+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(wind.ne, wind.se, ncol = 2)

### modis sst

load("e:/Data Manipulation/xtractomatic/sst.RData")
lon = sst.mod$longitude
lat = sst.mod$latitude
time = sst.mod$time %>% as.Date()
sst = sst.mod$data

dimension = data.frame(lon, sst[,,1]) %>% dim()
sst.kimbiji = list()

for (j in 1:length(time)){
  
sst.kimbiji[[j]] = data.frame(lon, sst[,,j]) %>% 
  gather(key = "lata", value = "sst", 2:dimension[2]) %>% 
  mutate(lat = rep(lat, each = dimension[1]), date = time[j], 
         month = month(date), season = month, 
         season = replace(season,season %in% c(11,12,1,2,3, 4), "NE"),
         # season = replace(season,season %in% c(4,10), "IN"),
         season = replace(season,season %in% c(5,6,7,8,9,10 ), "SE")) %>% 
  filter(lon >=38 & lon <= 42 & lat >=-8 & lat <= -3) %>%
  as.tibble()
}

## subset for the Area of interest
sst.kimbiji = sst.kimbiji %>% bind_rows() 

sst.kimbiji.monthly = sst.kimbiji %>% group_by(lon,lat, month) %>% summarise(sst = mean(sst, na.rm = TRUE))
sst.kimbiji.season = sst.kimbiji %>% group_by(lon,lat, season) %>% summarise(sst = mean(sst, na.rm = TRUE))
sst.modis.se = ggplot() +
  metR::geom_contour_fill(data = sst.kimbiji.season %>% filter(season == "SE" & sst <= 28), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_contour2(data = sst.kimbiji.season %>% filter(season == "SE" & sst <= 28), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE) + 
  metR::geom_text_contour(data = sst.kimbiji.season %>% filter(season == "SE" & sst <= 28), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, 
                          parse = TRUE, check_overlap = TRUE, size = 3) + 
  # metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
  #                   arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
  #                   pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColors9A(120), 
                       # limits = c(0,1.6),
                       # breaks =seq(0.2,1.6,.3),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1),
        axis.text.y = element_blank())+ 
  # scale_mag(max = 1.6, name = expression(speed~(ms^{-1})))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))



sst.modis.ne = ggplot() +
  metR::geom_contour_fill(data = sst.kimbiji.season %>% filter(season == "NE" & sst >28.1 & sst <= 30), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, bins = 70) + 
  metR::geom_contour2(data = sst.kimbiji.season %>% filter(season == "NE" & sst >28.1 & sst <= 30), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE) + 
  metR::geom_text_contour(data = sst.kimbiji.season %>% filter(season == "NE" & sst >28.1 & sst <= 30), 
                          aes(x = lon, y = lat, z = sst), na.fill = TRUE, 
                          parse = TRUE, check_overlap = TRUE, size = 3) + 
  # metR::geom_vector(data = uv.se, aes(x = lon, y = lat, dx = u, dy = v), 
  #                   arrow.angle = 20, arrow.type = "open", arrow.length = .5, 
  #                   pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Current",
                       colours = oceColors9A(120), 
                       # limits = c(0,1.6),
                       breaks =seq(28.2,30,0.4),
                       na.value = "white")+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  # scale_mag(max = 1.6, name = expression(speed~(ms^{-1})))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(sst.modis.ne, sst.modis.se, ncol = 2)

require(gganimate)

sst.kimbiji.monthly$month = sst.kimbiji.monthly$month %>% as.integer
wind.data.monthly$month = wind.data.monthly$month %>% as.integer


ggplot() +
  geom_raster(data = sst.kimbiji.monthly, 
                          aes(x = lon, y = lat, fill = sst), interpolate = TRUE) +
  # metR::geom_contour_fill(data = sst.kimbiji.monthly, 
  #                         aes(x = lon, y = lat, z = sst), na.fill = TRUE) + 
  metR::geom_vector(data = wind.data.monthly, 
                    aes(x = lon, y = lat, dx = wind.x/10, dy = wind.y/10), 
                    arrow.angle = 30, arrow.type = "open", arrow.length = .75, 
                    pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  scale_fill_gradientn(name = "Temperature",
                       colours = oceColorsJet(120), 
                       # limits = range(sst), 
                       # breaks =seq(0.2,1.6,.3), 
                       na.value = "white")+
  labs(x = "", y = "")+
 theme_bw()+
  theme(legend.position = "right",
        legend.key.height = unit(1.4, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_mag( name = expression(speed~(ms^{-1})))+
  labs(x = NULL, y = NULL, title = "Month of : {frame_time}")+
  transition_time(month) +
  ease_aes("linear")

ARGO float temperature and salinity

# list the netcdf file names in the working directory

argo.file = dir("e:/Doctoral/udsm/Processing/argo_profile/",  #relative path of the directory
                recursive = TRUE, # recurse into directory
                pattern = "_prof.nc", # choose the names match this condition
                full.names = TRUE) # get name with its relative path

## obtain id for files
argo.id = argo.file %>% 
  as.tibble() %>% 
  rename(profile = 1) %>% 
  separate(profile, c(NA, "id"), sep = "//") %>% 
  separate(id, c(NA, "Id", NA), sep = "/") %>% 
  pull(Id)
# preallocate a file to store processed ctd data
argo.ctd = NULL

# the first section of the loop run with i through the netcdf files
    for (i in 1:length(argo.file)){
        # read the files in sequence
      argo = read.argo(argo.file[i])%>%handleFlags()%>%argoGrid(p = seq(0,1000,5))
        # convert the argo list data into section
      argo.section = argo%>%as.section()
        # convert argo.section into list
      argo.list = argo.section[["station"]]

# the second section of the loop run with j through each argo list created by i
    for (j in 1:length(argo.list)){
      
    profile = argo.list[[j]]@data%>% # get each station profile data
      as.data.frame()%>%          # convert the data into daa frame
            # add variable of argo id, note this use i and not j
      mutate(id = argo.id[i], 
              #add station variable
             station = argo.list[[j]][["station"]],
              #add date of sampling variable
             time = argo.list[[j]][["startTime"]]%>%as.Date(),
              #add longitude of station variable
             lon = argo.list[[j]][["longitude"]], 
              #add latitude of station variable
             lat = argo.list[[j]][["latitude"]], 
              #compute the salinity variable in each station
             density = argo.list[[j]]%>%swRho(eos = "gsw"))%>%
        # convert the data frame into tibble
      as.tibble()%>%
      # drop other variable of no interest
      select(id, station, time, lon, lat,scan, pressure, 
             salinity = salinityAdjusted, temperature = temperatureAdjusted)
    #
    # # separate the wmoid into into diferent variables
    # profile = profile%>%separate(argo, c(1:7,"wmoid",8), sep = "/")%>%
    #   # keep the wmoid variable and drop the rest
    #   select(-c(1:4,6:10)) %>% rename(ID = 1)
    #
    #bind the argo.ctd with the data by rows
    argo.ctd = argo.ctd%>%bind_rows(profile)
  
}
}

ctd profiles

ctd collected during 2018 cruise

file = dir("e:/Data Manipulation/ctd_bongo2018/", pattern = ".cnv",full.names = TRUE)

ctd.se = list()

for (i in 1:length(file)){
  
  ctd = read.ctd(file[i])%>%ctdTrim(method = "downcast")%>%ctdDecimate(p = 5)
  
  ctd.se[[i]] = ctd@data %>% 
    as.tibble() %>% 
    mutate(time = ctd@metadata$startTime %>% as.Date(), 
           scan = i, station = i, id = ctd@metadata$station, 
           lon = ctd@metadata$longitude, lat = ctd@metadata$latitude)
  
}

ctd.se = ctd.se %>% bind_rows() %>% 
    select(id, station, time, lon, lat, scan, pressure, salinity, temperature)

ctd collected during the 2017 cruise

file = dir("e:/Data Manipulation/ctd_bongo2017/Converted Data/Bongo/", pattern = ".cnv",full.names = TRUE)


ctd.ne = list()

for (i in 1:length(file)){
  
  ctd = read.ctd(file[i])%>%
    ctdTrim(method = "downcast")%>%
    ctdDecimate(p = 5)
  
  ctd.ne[[i]] = ctd@data %>% 
    as.tibble() %>% 
    mutate(time = ctd@metadata$startTime %>% as.Date(), 
           scan = i, station = i, id = ctd@metadata$station, 
           lon = ctd@metadata$longitude, lat = ctd@metadata$latitude)
  
}

ctd.ne = ctd.ne %>% bind_rows() %>% 
    select(id, station, time, lon, lat, scan, pressure, salinity, temperature)
file = dir("e:/Data Manipulation/ctd_algoa/", pattern = "dstn",full.names = TRUE)


ctd.algoa = list()

for (i in 1:length(file)){
  
  ctd = read.ctd(file[i])%>%
    ctdTrim(method = "downcast")%>%
    ctdDecimate(p = 5)
  
  ctd.algoa[[i]] = ctd@data %>% 
    as.tibble() %>% 
    mutate(time = ctd@metadata$startTime %>% as.Date(), 
           scan = i, station = i, id = ctd@metadata$station, 
           lon = ctd@metadata$longitude, lat = ctd@metadata$latitude)
  
}

ctd.algoa = ctd.algoa %>% bind_rows() %>% 
    select(id, station, time, lon, lat, scan, pressure, salinity, temperature)
## because of the long time processing, i switched the chunk above not to process every time and instead load the products processed with them from previous session
load("argo_ctd1.RData")
## combine argo with ctd data from two Agulhas cruise and algoa


argo.ctd$station = as.integer(argo.ctd$station)
all.ctd = argo.ctd %>% bind_rows(ctd.ne, ctd.se, ctd.algoa)

## assign the seaoson
argo.ctd.season = all.ctd %>% mutate(month = month(time),
                    season = month, 
                    season = replace(season,season %in% c(11,12,1,2,3, 4), "NE"),
                    # season = replace(season,season %in% c(4,10), "IN"),
                    season = replace(season,season %in% c(5,6,7,8,9,10 ), "SE")) %>% na.omit()
## subset for the Area of interest

profile.map = argo.ctd.season %>% 
  filter(id ==  1900409 & pressure == 10 & lat > -9 & lon < 41) %>% 
  ggplot() + 
  geom_path(aes(lon,lat, group = id, col = season), size = .8) +
  geom_point(aes(lon,lat, col = season), size = 3) +
  
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-8.5,-5.5), xlim = c(38.5, 40.5))+
  ggrepel::geom_text_repel(aes(lon,lat,label = time%>%as.character()), size = 3.5)+
  theme_bw()+
  theme(legend.position = c(.2,.1),
        legend.key.height = unit(0.5, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 10, colour = 1),
        legend.text = element_text(size = 10),
        legend.key.width = unit(.85, "cm"),
        legend.key = element_blank(),
        legend.title = element_blank())+
  scale_x_continuous(breaks = c(38.7,40.2))+
  scale_y_continuous(breaks = c(-8.2, -5.7))+
  labs(x = NULL, y = NULL)


profile.temperature = argo.ctd.season %>% 
  filter(id ==  1900409 & pressure < 500 & lat > -9 & lon < 41) %>% 
  ggplot() + 
  geom_contour_fill(aes(x = lat, y = pressure, z = temperature), na.fill = TRUE, bins = 120)+
  geom_contour2(aes(x = lat, y = pressure, z = temperature), na.fill = TRUE)+
  geom_text_contour(aes(x = lat, y = pressure, z = temperature), 
                    na.fill = TRUE, check_overlap = TRUE, parse = TRUE, 
                    size =3.5, rotate = FALSE)+
  scale_y_reverse()+
 scale_x_continuous(breaks = seq(-8.2,-5.8, length.out = 4), 
                     labels = LatLabel(seq(-8.2,-5.8, length.out = 4)))+
  scale_fill_gradientn(colours = oce.colorsTemperature(120))+
  theme(legend.position = "none", axis.text = element_text(size = 10, colour = 1),
        panel.background = element_rect(colour = NA, fill = NA),
        axis.title = element_blank(), axis.text.y = element_blank(),
        axis.ticks.y = element_blank())+
  labs(x = NULL, y = NULL)+
  annotate(geom = "text", x = -7.8, y = 30, label = "Temperature")

profile.salinity = argo.ctd.season %>% 
  filter(id ==  1900409 & pressure < 500 & lat > -9 & lon < 41) %>% 
  ggplot() + 
  geom_contour_fill(aes(x = lat, y = pressure, z = salinity), na.fill = TRUE, bins = 120)+
  geom_contour2(aes(x = lat, y = pressure, z = salinity), na.fill = TRUE)+
  geom_text_contour(aes(x = lat, y = pressure, z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  scale_y_reverse(breaks = seq(50,450,100))+
  scale_x_continuous(breaks = seq(-8.2,-5.8, length.out = 4), 
                     labels = LatLabel(seq(-8.2,-5.8, length.out = 4)))+
  scale_fill_gradientn(colours = oce.colorsSalinity(120))+
  theme(legend.position = "none", axis.text = element_text(size = 10, colour = 1),
        panel.background = element_rect(colour = NA, fill = NA),
        axis.title.y = element_text(colour = 1, size = 12), plot.tag = element_text())+
  labs(x = NULL, y = "Water depth [meter]")+
  annotate(geom = "text", x = -8.1, y = 30, label = "Salinity")

egg::ggarrange(profile.map,profile.salinity, profile.temperature, nrow = 1)

creating polygons as AOI in sf from a series of lat and lon coordinates.

st_as_sf() has an argument coords that take points given as coordinate columns in a data frame and convert those columns to sf POINT geometries. Then, because sf works well with dplyr, we can use st_combine() to combine the points into a MULTIPOINT and use st_cast() to convert to POLYGON.

argo.ctd.season = argo.ctd.season %>% 
  filter(lon >= 37 & lon < 42 & lat > -7 & lat < -4)

argo.ctd.sf = argo.ctd.season %>% 
  st_as_sf(coords = c("lon", "lat")) %>% 
  st_set_crs(4326)



# extent.df = data.frame(lon =c(38.47412,38.47412, 50.44922, 50.44922,38.47412), 
#                        lat = c(-13.20576,0.350621, 0.350621, -13.20576,-13.20576))
# 
# extent= extent.df %>% 
#   st_as_sf(coords = c("lon", "lat")) %>%
#   st_set_crs(4326) %>% 
#   summarise(geometry = st_combine(geometry)) %>%
#   st_cast("POLYGON")
grid = argo.ctd.sf %>% st_make_grid(n = 20) %>% st_sf()

ggplot() + 
  geom_sf(data = grid)+
  geom_sf(data = tz.ke, fill = "grey95", col = 1)+
  coord_sf(xlim = argo.ctd.season$lon %>% range, ylim = argo.ctd.season$lat %>% range)

northeast season

argo.ctd.sf.ne = argo.ctd.sf %>% filter(pressure == 10 & season == "NE") 

argo.grid.ne= grid %>%
  mutate(id = 1:n(),
         contained = lapply(st_contains(st_sf(geometry), argo.ctd.sf.ne), identity),
         obs = sapply(contained, length), 
         temperature = sapply(contained, function(x) {mean(argo.ctd.sf.ne[x,]$temperature, na.rm = TRUE)}),
         salinity = sapply(contained, function(x) {mean(argo.ctd.sf.ne[x,]$salinity, na.rm = TRUE)}))

argo.grid.ne = argo.grid.ne %>% mutate(pressure = 10)%>%select(pressure,obs, temperature,salinity)

depth = c(100,200,500,1000)

for (i in depth){
  
  argo.ctd.sf.ne = argo.ctd.sf %>% filter(pressure == i & season == "NE") 

argo.sf.ne = grid %>%
  mutate(id = 1:n(),
         contained = lapply(st_contains(st_sf(geometry), argo.ctd.sf.ne), identity),
         obs = sapply(contained, length), 
         temperature = sapply(contained, function(x) {mean(argo.ctd.sf.ne[x,]$temperature, na.rm = TRUE)}),
         salinity = sapply(contained, function(x) {mean(argo.ctd.sf.ne[x,]$salinity, na.rm = TRUE)}))

argo.sf.ne = argo.sf.ne %>% mutate(pressure = i)%>%select(pressure,obs, temperature,salinity)
  
argo.grid.ne = argo.grid.ne %>% rbind(argo.sf.ne)
}
argo.ne.df = argo.grid.ne
coordinates = argo.grid.ne %>% st_centroid() %>% st_coordinates() %>% as.tibble() %>% rename(lon=1,lat=2)
st_geometry(argo.ne.df) = NULL

argo.ne.df = coordinates %>% bind_cols(argo.ne.df) %>% na.omit()

temprature

depth = c(10, 100, 200, 500, 1000)


ne.df.temperature = NULL

for (n in depth) {
  
ten = argo.ne.df %>% filter(pressure == n)

xg = pretty(x = ten$lon, n = 40)
yg = pretty(x = ten$lat, n = 40)

 argo.ne.interp = interpBarnes(ten$lon, ten$lat, ten$temperature, xg = xg, yg = yg)
 
 dim.lon = length(argo.ne.interp$xg)
dim.lat = length(argo.ne.interp$yg)

 argo.ne.interp.df = data.frame(argo.ne.interp$xg,  argo.ne.interp$zg%>%as.data.frame())%>%
   as.tibble()%>%
   gather(key = "aa", value = "temperature", 2:(dim.lat+1))%>%
   mutate(lat = rep(argo.ne.interp$yg, each = dim.lon), pressure = n)%>%
   select(pressure, lon = 1, lat, temperature)%>%na.omit()
 
ne.df.temperature = ne.df.temperature %>% bind_rows(argo.ne.interp.df) 
}

salinity

depth = c(10, 100, 200, 500, 1000)


ne.df.salinity = NULL

for (n in depth) {
  
ten = argo.ne.df %>% filter(pressure == n)

xg = pretty(x = ten$lon, n = 40)
yg = pretty(x = ten$lat, n = 40)

 argo.ne.interp = interpBarnes(ten$lon, ten$lat, ten$salinity, xg = xg, yg = yg)
 
 dim.lon = length(argo.ne.interp$xg)
dim.lat = length(argo.ne.interp$yg)

 argo.ne.interp.df = data.frame(argo.ne.interp$xg,  argo.ne.interp$zg%>%as.data.frame())%>%
   as.tibble()%>%
   gather(key = "aa", value = "salinity", 2:(dim.lat+1))%>%
   mutate(lat = rep(argo.ne.interp$yg, each = dim.lon), pressure = n)%>%
   select(pressure, lon = 1, lat, salinity)%>%na.omit()
 
ne.df.salinity = ne.df.salinity %>% bind_rows(argo.ne.interp.df) 
}

ne.df.all = ne.df.salinity %>% bind_cols(ne.df.temperature %>% select(temperature))
d0.temperature = ggplot()+
  geom_contour_fill(data = ne.df.all %>% filter(pressure == 10),
                    aes(x = lon, y = lat, z = temperature), na.fill = TRUE, bins = 70)+
  geom_contour2(data = ne.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE)+
  geom_text_contour(data = ne.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
 coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

d200.temperature = ggplot()+
  geom_contour_fill(data = ne.df.all %>% filter(pressure == 200),
                    aes(x = lon, y = lat, z = temperature), na.fill = TRUE, bins = 70)+
  geom_contour2(data = ne.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE)+
  geom_text_contour(data = ne.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(d0.temperature, d200.temperature, ncol = 2)

d0.salinity = ggplot()+
  geom_contour_fill(data = ne.df.all %>% filter(pressure == 10),
                    aes(x = lon, y = lat, z = salinity), na.fill = TRUE, bins = 70)+
  geom_contour2(data = ne.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE)+
  geom_text_contour(data = ne.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

d200.salinity = ggplot()+
  geom_contour_fill(data = ne.df.all %>% filter(pressure == 200),
                    aes(x = lon, y = lat, z = salinity), na.fill = TRUE, bins = 70)+
  geom_contour2(data = ne.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE)+
  geom_text_contour(data = ne.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+ 
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1),
        axis.text.y = element_blank())+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(d0.salinity, d200.salinity, ncol = 2)

southeast monsoon

argo.ctd.sf.se = argo.ctd.sf %>% filter(pressure == 10 & season == "SE") 

argo.grid.se= grid %>%
  mutate(id = 1:n(),
         contained = lapply(st_contains(st_sf(geometry), argo.ctd.sf.se), identity),
         obs = sapply(contained, length), 
         temperature = sapply(contained, function(x) {mean(argo.ctd.sf.se[x,]$temperature, na.rm = TRUE)}),
         salinity = sapply(contained, function(x) {mean(argo.ctd.sf.se[x,]$salinity, na.rm = TRUE)}))

argo.grid.se = argo.grid.se %>% mutate(pressure = 10)%>%select(pressure,obs, temperature,salinity)

depth = c(100,200,500,1000)

for (i in depth){
  
  argo.ctd.sf.se = argo.ctd.sf %>% filter(pressure == i & season == "SE") 

argo.sf.se = grid %>%
  mutate(id = 1:n(),
         contained = lapply(st_contains(st_sf(geometry), argo.ctd.sf.se), identity),
         obs = sapply(contained, length), 
         temperature = sapply(contained, function(x) {mean(argo.ctd.sf.se[x,]$temperature, na.rm = TRUE)}),
         salinity = sapply(contained, function(x) {mean(argo.ctd.sf.se[x,]$salinity, na.rm = TRUE)}))

argo.sf.se = argo.sf.se %>% mutate(pressure = i)%>%select(pressure,obs, temperature,salinity)
  
argo.grid.se = argo.grid.se %>% rbind(argo.sf.se)
}
argo.se.df = argo.grid.se
coordinates = argo.grid.se %>% st_centroid() %>% st_coordinates() %>% as.tibble() %>% rename(lon=1,lat=2)
st_geometry(argo.se.df) = NULL

argo.se.df = coordinates %>% bind_cols(argo.se.df) %>% na.omit()

temprature

depth = c(10, 100, 200, 500, 1000)


se.df.temperature = NULL

for (n in depth) {
  
ten = argo.se.df %>% filter(pressure == n)

xg = pretty(x = ten$lon, n = 40)
yg = pretty(x = ten$lat, n = 40)

 argo.se.interp = interpBarnes(ten$lon, ten$lat, ten$temperature, xg = xg, yg = yg)
 
 dim.lon = length(argo.se.interp$xg)
dim.lat = length(argo.se.interp$yg)

 argo.se.interp.df = data.frame(argo.se.interp$xg,  argo.se.interp$zg%>%as.data.frame())%>%
   as.tibble()%>%
   gather(key = "aa", value = "temperature", 2:(dim.lat+1))%>%
   mutate(lat = rep(argo.se.interp$yg, each = dim.lon), pressure = n)%>%
   select(pressure, lon = 1, lat, temperature)%>%na.omit()
 
se.df.temperature = se.df.temperature %>% bind_rows(argo.se.interp.df) 
}

salinity

depth = c(10, 100, 200, 500, 1000)


se.df.salinity = NULL

for (n in depth) {
  
ten = argo.se.df %>% filter(pressure == n)

xg = pretty(x = ten$lon, n = 40)
yg = pretty(x = ten$lat, n = 40)

 argo.se.interp = interpBarnes(ten$lon, ten$lat, ten$salinity, xg = xg, yg = yg)
 
 dim.lon = length(argo.se.interp$xg)
dim.lat = length(argo.se.interp$yg)

 argo.se.interp.df = data.frame(argo.se.interp$xg,  argo.se.interp$zg%>%as.data.frame())%>%
   as.tibble()%>%
   gather(key = "aa", value = "salinity", 2:(dim.lat+1))%>%
   mutate(lat = rep(argo.se.interp$yg, each = dim.lon), pressure = n)%>%
   select(pressure, lon = 1, lat, salinity)%>%na.omit()
 
se.df.salinity = se.df.salinity %>% bind_rows(argo.se.interp.df) 
}

se.df.all = se.df.salinity %>% bind_cols(se.df.temperature %>% select(temperature))
d0.temperature = ggplot()+
  geom_contour_fill(data = se.df.all %>% filter(pressure == 10),
                    aes(x = lon, y = lat, z = temperature), na.fill = TRUE, bins = 70)+
  geom_contour2(data = se.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE)+
  geom_text_contour(data = se.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  # metR::geom_vector(data = uv.se %>% sample_frac(.25), aes(x = lon, y = lat, dx = u, dy = v), 
  #                   arrow.angle = 40, arrow.type = "open", arrow.length = .75, 
  #                   pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  # scale_mag(max = 1.8, name = expression(speed~(ms^{-1})))+
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

d200.temperature = ggplot()+
  geom_contour_fill(data = se.df.all %>% filter(pressure == 200),
                    aes(x = lon, y = lat, z = temperature), na.fill = TRUE, bins = 70)+
  geom_contour2(data = se.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE)+
  geom_text_contour(data = se.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = temperature), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
  # metR::geom_vector(data = uv.se %>% sample_frac(.5), aes(x = lon, y = lat, dx = u, dy = v), 
  #                   arrow.angle = 40, arrow.type = "open", arrow.length = .75, 
  #                   pivot = 0,preserve.dir = TRUE, direction = "ccw")+
  geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"),
        legend.background = element_blank(),
        axis.text = element_text(size = 12, colour = 1))+ 
  # scale_mag(max = 1.8, name = expression(speed~(ms^{-1})))+
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~degree~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(d0.temperature, d200.temperature, ncol = 2)

d0.salinity = ggplot()+
  geom_contour_fill(data = se.df.all %>% filter(pressure == 10),
                    aes(x = lon, y = lat, z = salinity), na.fill = TRUE, bins = 70)+
  geom_contour2(data = se.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE)+
  geom_text_contour(data = se.df.all %>% filter(pressure ==10), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
     geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"), 
        legend.background = element_blank(),
        axis.text = element_text(size = 11, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~Salinity~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-5.9, -4.4))

d200.salinity = ggplot()+
  geom_contour_fill(data = se.df.all %>% filter(pressure == 200),
                    aes(x = lon, y = lat, z = salinity), na.fill = TRUE, bins = 70)+
  geom_contour2(data = se.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE)+
  geom_text_contour(data = se.df.all %>% filter(pressure ==200), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
   geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"),
        legend.background = element_blank(), axis.text.y = element_blank(),
        axis.text = element_text(size = 11, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~Salinity~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-5.9, -4.4))

d500.salinity = ggplot()+
  geom_contour_fill(data = se.df.all %>% filter(pressure == 500),
                    aes(x = lon, y = lat, z = salinity), na.fill = TRUE, bins = 70)+
  geom_contour2(data = se.df.all %>% filter(pressure ==500), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE)+
  geom_text_contour(data = se.df.all %>% filter(pressure ==500), 
                aes(x = lon, y = lat , z = salinity), na.fill = TRUE, 
                    check_overlap = TRUE, parse = TRUE, size =3.5, rotate = FALSE)+
   geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(1.2, "cm"), 
        legend.key.height = unit(0.2, "cm"),
        legend.background = element_blank(), axis.text.y = element_blank(),
        axis.text = element_text(size = 11, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~Salinity~C))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-5.9, -4.4))

egg::ggarrange(d0.salinity, d200.salinity, ncol = 2)

AVISO-SLA

variables = metR::ReadNetCDF("e:/MatlabWorking/Altimetry/msla_h/indian_ocean-twosat-msla-h_010193_311295.nc", out = "data.frame")

 variables$dimensions %>% names()
 
sla = metR::ReadNetCDF("e:/MatlabWorking/Altimetry/msla_h/indian_ocean-twosat-msla-h_010193_311295.nc", subset = list(lon= -35:43, level = as.Date("1995-06-10")))
nc = ncdf4::nc_open("e:/MatlabWorking/Altimetry/msla_h/indian_ocean-twosat-msla-h_010193_311295.nc")
time = ncdf4::ncvar_get(nc, "time")
lon = ncdf4::ncvar_get(nc, "lon")
lat = ncdf4::ncvar_get(nc, "lat")
sla  = ncdf4::ncvar_get(nc, "sla")

sla %>% dim; lon %>% length(); lat %>% length()
to = ymd("1950-01-01", tz = "")
to.jd = insol::JD(to)
t = time+to.jd
t = insol::JD(t, inverse = TRUE) 
sla.list = list()

week.day = seq(1,length(t),8)

dimension = data.frame(lon, sla[,,1]) %>% dim()

# for (j in 1:length(t)){
for (j in week.day){
sla.list[[j]] = data.frame(lon, sla[,,j]) %>% 
  gather(key = "lata", value = "sla", 2:dimension[2]) %>% 
  mutate(lat = rep(lat, each = dimension[1]), time = t[j]) %>% 
  separate(time, c("time", NA), sep = " ") %>% 
  filter(lon > 37 & lon < 45 & lat > -14 & lat < -1)
}

sla.tb = sla.list %>% bind_rows() 
load("sla.RData")
sla.tb.season = sla.tb %>%
  mutate(month = month(time),
         season = month, 
         season = replace(season,season %in% c(11,12,1,2,3, 4), "NE"),
         # season = replace(season,season %in% c(4,10), "IN"),
         season = replace(season,season %in% c(5,6,7,8,9,10 ), "SE"),
         sla = sla*10) 
msla = sla.tb.season %>% group_by(lon,lat, month) %>% summarise(msla = mean(sla, na.rm = TRUE))
msla.pemba = msla %>% filter(lon >38.2 & lon < 41.2 & lat > -6.7 & lat < -4)

sla.ne = ggplot()+
  metR::geom_contour_fill(data = msla.pemba %>% filter(month==1),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE, bins = 120)+
  metR::geom_contour2(data = msla.pemba %>% filter(month==1),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE)+
  metR::geom_text_contour(data = msla.pemba %>% filter(month==1),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE, 
                    check_overlap = FALSE, parse = TRUE, size =3.5, rotate = FALSE)+
   geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(.4, "cm"), 
        legend.key.height = unit(1.2, "cm"),
        legend.background = element_blank(), 
        axis.text = element_text(size = 11, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~SLA~(mm)))+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

sla.se = ggplot()+
  metR::geom_contour_fill(data = msla.pemba %>% filter(month==8),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE, bins = 120)+
  metR::geom_contour2(data = msla.pemba %>% filter(month==8),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE)+
  metR::geom_text_contour(data = msla.pemba %>% filter(month==8),
       aes(x = lon, y = lat, z = msla), na.fill = TRUE, 
                    check_overlap = FALSE, parse = TRUE, size =3.5, rotate = FALSE)+
   geom_sf(data = tz.ke,fill = "lightgrey", col = "black")+
  coord_sf(ylim = c(-6.5,-4.5), xlim = c(38.5, 41))+
  labs(x = "", y = "")+
  theme_bw()+
  theme(legend.position = "none",
        legend.key.width = unit(.4, "cm"), 
        legend.key.height = unit(1.2, "cm"),
        legend.background = element_blank(), axis.text.y = element_blank(),
        axis.text = element_text(size = 11, colour = 1))+ 
  scale_fill_gradientn(colours = oceColorsJet(120), name = expression(~SLA~(mm)))+
  # scale_fill_divergent(midpoint = 0)+
  scale_x_continuous(breaks = c(38.8,40.5))+
  scale_y_continuous(breaks = c(-6.3, -4.8))

egg::ggarrange(sla.ne, sla.se, ncol = 2)
Mean Sea Level anomaly during a) NE and b) SE monsoon season

Mean Sea Level anomaly during a) NE and b) SE monsoon season